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The inspiral and merger of supermassive black hole binary systems with high orbital eccentricity 
are among the promising sources of the advanced gravitational wave observatories. In this paper 
we derive analytic ready-to-use first post-Newtonian eccentric waveform in Fourier domain with 
the use of Hansen coefficients. Introducing generic perturbations of celestial mechanics we have 
generalized the Hansen expansion to the hrst post-Newtonian order which are then used to express 
the waveforms. Taking into account the high eccentricity of the orbit leads to the appearance of 
secular terms in the waveform which are eliminated with the introduction of a phase shift. The 
waveforms have a systematic structure and as our main result these are expressed in a tabular form. 
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I. INTRODUCTION 

Compact binaries (i. e. black holes, neutron stars and white dwarfs) with non-vanishing eccentricity are promising 
sources of gravitational waves. Depending on their parameters these sources emit radiation within the sensitivity 
band of the forthcoming gravitational wave detectors advanced LIGO and Virgo. Due to their increased sensitivity 
the source signals will be visible by these detectors for a longer time period with the requirement of an accurate 
description of both the orbital evolution and the emitted GWs of these systems, detectability The measured signal 
output of the detectors is cross correlated with theoretical waveform templates in matched filtering. The presence of 
orbital eccentricity changes significantly the properties of the waveforms resulting in the decrease of their detectability 
with the use of circular templates.In spite of the general circularization of binary orbits due to GW emission binaries 
interacting with their environment can retain non-negligible eccentricity at the end of their evolution. For example, 
there are indications that binaries in dense galactic nuclei [H, , embedded in a gaseous disk ii can remain eccentric 

until the end of their inspiral. Moreover, the interaction of supermassive black hole binaries with star populations 
[^, @ and the Kozai mechanism and relativistic orbital resonances in hierarchial triples [M3 can also increase orbital 
eccentricity. 

The first analytic eccentric waveform for the unperturbed motion was given by m- The time-dependent waveform 
is computed with the help of Fourier-Bessel method in Ref. [l^ and the waveform in Fourier space for arbitrary 
eccentricity was given in Ref (T^ . 

Recently, within the post-Newtonian (PN) treatment of compact binary evolution the theoretical computations are 
reachiM the level of 4PN. The first post-Newtonian eccentric waveform for bound orbits is computed by Wagoner and 
Will [l4ll i n 1976. The description of the Kepler motion with the IPN correction is given by Damour and Deruelle 
in Ref [l5l| which requires three eccentricities {radial, time and angular eccentricity). With the help of this Damour- 
Deruelle parameterization the eccentric waveform and evolution of the semimajor axis and the radial eccentricity due 
to radiation reaction is computed by Junker and Schafer in Ref [T^ . The waveform in Fourier domain up to IPN and 
2PN orders are given in Refs. [13,H- 

Our work focuses on ready-to-use eccentric IPN waveforms in frequency domain. In our description we use the 
generalized true anomaly parameterization which has the advantage that the solution of the equations of motion can 
be expressed with only two eccentricities. As a consequence, the gravitational waveforms have a simple structure. 
Secular terms appearing in the waveforms are eliminated by the use of the Poincare-Lindstedt method and the 
introduction of the drift true anomaly parameter. We give both the time-dependent and frequency domain waveforms 
using the Hansen expansion and the stationary phase approximation (SPA). Moreover, we compute the evolution of 
the semimajor axis and radial eccentricity up to IPN accuracy. 

The paper is organized as follows. After a short summary of the Damour-Deruelle parameterization in Sec. II 
we introduce the generalized true anomaly parameterization in Sec. Ill and the IPN waveform in Sec. IV. Sec. V. 
contains the extension of the Hansen coefficients to IPN order. Fourier domain SPA waveforms are given in Sec VI. 
The radiation reaction problem and the evolution of the time and phase functions to IPN order are given in Sec VH. 
Some of the technical details, i.e. tensor spherical harmonics, orbital parameters to IPN order, Hansen coefficients 
and waveform expressions are presented in the Appendices. 


II. DAMOUR-DERUELLE PARAMETERIZATION 


In the following, we summarize the first post-Newtonian parameterization of the orbital motion introduced by 
Damour and Deruelle [l^] for the description of compact binaries. The Lagrangian of the system contains the 
Newtonian and IPN corrections, 
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where r is the relative distance, v is the relative velocity vector, m = mi -I- m 2 is the total mass, /x = mim 2 /m is the 
reduced mass and rj = /i/m is the symmetric mass ratio. The overdot denotes the derivative with respect to time t. 
The radial and angular motion can be separated to the linear order, so the Euler-Lagrange equation are 
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where 9 is the azimuthal angle in the orbital plane and the constants A, B,C, D, H, I depend on the conserved 
quantities such as energy and the orbital angular momentum of the perturbed motion, see Appendix B. The constants 
A, B, C, H contain Newtonian and IPN terms while D and I are purely IPN corrections. The equations of motion, 
Eqs. (I2I3I) can be solved by the eccentric anomaly quasi-parameterization u, that is 

r = 0^(1 — Cj. cosu) , (4) 

where the orbital parameters are the semimajor axis Or and the radial eccentricity e^. These orbital parameters are 
characterized by the turning points (rmax and rmin in [2^,[l^) of the radial motion. The Kepler equation and angular 
evolution can be given as 


nit — to) = u — et sin u = A4, 

(5) 

9 - 9o = {1 + k)vg , 

( 6 ) 
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in terms of the orbital elements of the IPN orbital dynamics such as the mean motion n, the time eccentricity et, the 
angle eccentricity eg, and the pericenter drift k (which is in relationship with the pericenter precession ( 7 ) averaged 
over one radial period, see [1^). In the above equations to and 9o are integration constants and in our calculations 
we set 9o = to = 0. The Damour-Deruelle parameterization contains three eccentricities, which are different in IPN 
order only. In Newtonian order only the Kepler equation and one eccentricity remain. 

We note that the Kepler equation remains the same with the inclusion of the spin-orbit interaction (see ref. [Hi) 
but it contains higher order contributions such as the spin-spin, quadrupole-monopole and magnetic dipole-dipole 


interactions [20j or the second post-Newtonian corrections [22|. For this reason we use an other parameterization 
which has only two eccentricities and the evolution of the azimuthal angle 9 is similarly simple as Eq. ®. 


III. GENERALIZED TRUE ANOMALY PARAMETERIZATION 


We introduce the generalized true anomaly parameterization (denoted by y in 


[2l,[23) 


as 


arjl - el) 
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This parameterization has a similar form than the Keplerian parameterization, where the orbital parameters 
Qr and er include only the leading order Newtonian terms. Eqs. ([ 8 ]) and o lead to the relation (j) = Vr = 

2arctan (1-|-(1 — tan (u/2) , cf. Eq. ([71). 

The time evolution of the generalized true anomaly up to IPN order can be expressed as 

na^(l - 

r'^il — etCr + (cr — et) cos tf) 

The integration of Eq. ([3]) with the help of (0 leads to the relation 

9 — 9o = (1 + Ki) + K 2 sin 0, 


(9) 


( 10 ) 


where we have introduced the IPN order quantities Ki = 3Gm/[af.(l — e^)c^] and K 2 = Gp,er/[2ar{l — e^)c^]. In this 
parameterization only radial and time eccentricities {Cr-, et) appear while angle eccentricity does not. The relationship 
between (j) and vg is given by Eqs. m and (HOD up to IPN order as 


(j) = vg - 


Gp,er 

2ar(l — e^)c^ 


sin Vg. 


( 11 ) 


In the following we will compute the full eccentric IPN waveform with the use of the generalized true anomaly 
parameterization (j) without the appearance of secular terms in the expressions and give the time-domain waveforms 
using the generalized Hansen expansion. Our aim is to express the full analytic eccentric frequency-domain waveform 
up to IPN order. 
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IV. WAVEFORM 

The first explicit description of the IPN waveform for binaries was given by Wagoner and Will in 1976 [ij. We 
rewrite their expressions with help of Ref. [s^ and, as a result, the radiation field up to IPN order has the following 
form 
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Here is the luminosity distance, e = Gm/c^r is the post-Newtonian parameter, rj^^’^"*and are the tensorial 

electric and magnetic scalar harmonics which are given by Eq. (2.30d) in [30j. The quantities I km, S km are the 
kth time derivatives of the mass and current multipole moments. The explicit form of these multipoles was given 
by Junker and Schafer in IPN order with eccentric anomaly u in refs, [l^ and Later, the authors of [l^ have 

^ (fc) (fc) 

computed the explicit time-dependent multipoles I km and S km up to 2PN. Their waveforms contain no secular 
terms because the authors are not using Eq. (jS]) but the exponents containing 6 have been expressed as series in 

(fe) (fc) 

I km, S km ■ 

It is well-known that some secular terms will appear in the eccentric waveforms if one expands the harmonic 
functions of the angle 9 in terms of the generalized true anomaly parameterization. These secular terms have to be 
eliminated in the waveforms which requires the introduction of the drift true anomaly (j)' = {\ -\- ki) (j). Then the 
harmonic functions of 6 can be described in a perturbative sense (see [^) as 

cos 9 ~ cos (j) — K 2 sin (j) sin (p', 
sin 9 ~ sin (j)' + K 2 sin (f> cos (j)' > 

which relations will be used to eliminate secular terms in the eccentric waveforms. 

The waveform up to IPN order for the polarization states are 

h+,x{(t>) = ^?.x(</>) + ^+,x(</') + ( 16 ) 
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The X {(p) is the Newtonian, ^ (^) i® h®-!! order PN and {(p) is the IPN waveform (see Appendix D and 
E). Our IPN waveform has the well-known structure with the generalized true anomaly (p and drift anomaly (p', as 
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where the coefficients c'^^^y ®m(c 4 )’Cm(c 2 ) > '®m(c 4 )’ '^m^ ^ud depend on the radial eccentricity e^, the mass 
parameters m, fi, rj and the two polar angles 0 and of the line of sight (see Appendix E). We introduce the 
quantities 
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with the real numbers mi± = to ± 4 (1 + ki), TO 2 ± = to ± 2 (1 + ki) and coefficients 
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We can use the Hansen expansion for the generalized true anomalies as 

oo oo 

i A(^ = cos fcA4, sinA(^ = 5^ sin kM 
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with Cq = Xq’^, and where X]^'^ are the generalized Hansen coefficients for the 

real number A (since ki is not integer) up to IPN order. We note that in the Keplerian case (where A is integer and 
e = Br = Bt = eg) cosAf/) and sinA^ can be extended by trigonometric functions of cos(j) and sint/f. Using the Fourier 
coefficients containing Bessel-functions (see e.g. 
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where / denotes the derivative with respect to the eccentricity e. This ’classical’ extension leads to an increasing order 
of sums for the increasing value of A. Note that Eqs. ([Ml) and (113 are not valid for the IPN motion. Therefore we 
extend Hansen coefficients up to IPN order in the next chapter. 


V. GENERALIZATION OF HANSEN COEFFICIENTS 

The Hansen coefficients are well-known already since the 19th century in celestial mechanics (see Appendix C). In 
our description of the time-dependent waveforms there appear Hansen coefficients therefore it is important to extend 
the Hansen expansion up to IPN order. The Hansen coefficients appear in next series 
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The definition of Hansen coefficients are 
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In the waveform there appear the harmonic functions of toU where to is not an integer parameter. So we have to 
generalize the formula of the KBplBrian Hansen coefficients 27| , 
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It is evident if p is an integer {i.e. m is integer) then Ip{z) = Jp{z) where Jp{z) is the Bessel function {e.g. for the 
Newtonian waveform see [l3l|l. If p is not an integer then Ip{z) = Jp{z) + gp{z), where gp{z) is the correction integral 

(P) 


OO 
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for R{z) > 0. 

Due to the different eccentricities we have to generalize Hansen coefficient in a different way for case of IPN order. 
The mean anomaly (see Eq. (IC3I) in Appendix C) is 
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where Pr = P{^r) (see Appendix C, /3 = (1 — y/1 — e^)/e is the parameter in celestial mechanics). Then the integrand 
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which can be extended in a form of an infinite series of the sum. Then the generalized Hansen coefficients for IPN 
order are 
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If we use a new index notation j = t — s then generalized Hansen coefficients for IPN are 
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where p = k — m -\- j. We note that for an integer m the square bracket in second line of Eq. (1361) can be written as 
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Then the explicit time-dependent waveforms, Eq. dm) , are 
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These waveform has significantly simpler structure than the corresponding expressions in 0. 
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VI. WAVEFORM IN FOURIER SPACE 


The waveform in Fourier space can be described in the stationary phase approximation of the time-dependent 
waveform (see Eqs. (B2) and (B3) in the Appendix B of [l^)- Taking an arbitrary harmonic function A{t) cos<i>(t), 
where A(t), $(t) is the time-dependent amplitude and phase, respectively, and the conditions A/A <C $ and $ <C 
are satisfied, then the Fourier transform of the function A{t) cos$(t) can be written as 
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where 'f' [/(T)] = 27r/(T)t [j^(T)] — $ [i^{T)] is the phasing function, T is the saddle point and the functions t [j^(T)] 
and $ [*^(7^)] appearing in the above expressions can be obtained from the leading order equations for gravitational 
radiation by Appell functions (see the Appendix in El)- It is necessary to add, that here the phase and frequency ('I'n 
and 'I'„±) are not split ting into triplet due to pericenter precession (it was a consequence of the appearance of 7 , i.e. 
heuristic precession in |13|1. because it is contained directly the IPN equations of motion (see the orbital parameter 
k). Therefore the IPN waveform depends on the single phase and frequency Accordingly, the waveform, Eq. 
(l38l) . in the Fourier space becomes 
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with the stationary phase condition 
above form can be written as 
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where the phasing functions are = 27r/t(T) — 4>„(T) ± 7 r/ 4 . 

Afterwards we shall compute the phase <I>„(T) and time t(T) functions appearing in the IPN waveform. 


VII. RADIATION REACTION TO IPN ORDER 


To leading order the averaged radiative change of the semimajor axis and eccentricity is governed by the quadrupole 
formula, see Peters [ 2 ^. In these equations the semimajor axis can be replace by the orbital frequency using Kepler’s 
third law to have the following expressions 
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where A4c = mrf/^ is the chirp mass of the binary system. The above equations can be integrated and with the 
use of the exact solution the phase and time functions can be expressed in terms of the Appell functions (see [II). 
Afterwards, we will compute IPN corrections to these equations. 

The averaged losses of the radial orbital parameters due the gravitational radiation reaction up to IPN order is 
given by Junker and Schafer [T^. We have to use Kepler’s third law in IPN order relating the orbital frequency and 
semimajor axis as 


ar = A I 1 + (r; - 9) ^ 
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The radiative evolution of the orbital frequency and eccentricity up to IPN order can be written as (in this chapter 
we omit the subscript r of the radial eccentricity) 
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Thereafter we will find the perturbative solution to the above equations up to IPN order. 
We get the relation between p and e from Eqs. (l50l) and (ISTIl up to IPN order as 
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The exact general solution in the Newtonian order (without the two last terms in the right hand side of Eq. (IM1) 1 is 
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where pn is the general Newtonian solution. Hereafter we use the expression Pf^ = PQa{e)la{eo) where the quantities 

pq and 60 are initial values for PN{eo) = and a{e) = (l — (1 + ^ 6 ^)“^^°^/^^®® is a shorthand notation. 

The perturbative Eq. (IMl) has the exact general solution (including the Newtonian and IPN terms) in a mathematical 
sense 
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with the coefficients 

1153 89 

3192 ~ 

2293125927 , 60619 

~ 558758080 6984476’^’ 

86928802699 , 129501097^ 

”93871357440 670509696^’ 

703 785 517 49 913 

4014 235 680 ” 735 208’^' ^ ^ 

Here 2 ^i is the ordinary hypergeometric function. These general solutions for iy{e) and a(e) are consistent 

with the IPN order Kepler equation of (l4^ . strictly speaking if we solve similarly method the original radiation 
reaction (equations of {d, e} depending on the semimajor axis) equation up to IPN order with the use of Eq. (|4^ 
then we will get an identity for semimajor axis a (note that a = for Newtonian order). 

Let us identify the Newtonian expression = vjq = C'ocr(e), where Cq = vo/aieef). The integration constant 

C has leading order corrections at IPN order so if we require the equation v{eo) = vq to hold, we get the valid 
perturbative solution for the orbital frequency, Eq. (15611 . to IPN accuracy 
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up to IPN order. The integrals in the Newtonian case is given in Appendix A of 

15c^Ar 

304(GA4c) 5/3 27ri/o J 
re g/29/19(^_g/2)-3/2^g/ 

Jo (l+m^,2)-118V2299 


(63) 


($-$e)iv = - 


^ 304(GA4c)®/3 

l-e 

X 


-de , 


Jo (l + i§e'2)124/2299“"’ (^4) 

where we have introduced the notation Aq = CT(eo)/(27rr'o) which depends on the initial eccentricity and orbital 
frequency. Such type of integrals can be given by extended hypergeometric functions {i.e. Appell functions). 
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where Fi (a, /3,/3', 7 ; x, y) is the Appell function (see [H) and the constants are a = —10/19, d = 8/19, /3 = 
124/2299, ,5 = —1181/2299,7 = 3/2, <5 = —121/304 [/. Similiar integrands appear in IPN order. Then we can 


^ We have used the following integral formula for the Appell function 
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compute the integrand of time function to IPN order as 
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Note that in Eq. (IS71) the last term is coming from the solution of the orbital frequency in Eq. (|5(I1) . Compntation 
of the integral Eq. (1571) is difficult thus we use the approximation 2 E 1 - 5 ^, ff) ~ 1 because its limit is 1 

for e —>■ 1 and 0.9474 for e —>■ 1. Then the time function is 
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The final result is 
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where we have introduced the shorthand notations 7 = 3/2,5 = —121/304, /3 = 1988/2299, /3 = —1181/2299, 
tti = 1 — (ao + 2i), ao = 17/19 and Tq = Po^^bppf{eo). The phase function $ — $c can be computed similarly. The 
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integrand of the phase function $ — $c up to linear order is 

5c^Ao 


(‘f’ — 




76GMcV^/^ Jo (1 + i§ie'2)124/2299 

Bi + + Bse'^ 

(1 + i^e'2)3169/2299 


X 

(1 



+ 

27^1 

15^ 


304' 

870 13 32 121 

2299’19’19’”^® 


JN ■ 


/2 


- 646 ' 


/2 


de' 


(71) 


where we have introduced the quantities Bi = Bi — 9Bi/8. We use the above approximation 2 ^i(--) — !> thus the 
final form of the phase function is 


(*& ^c)p7v 


Sc^Aq 

76GM~n^ 








a2 



2 + ci;2 
2 



B 4 


ISFo 

2 




(72) 


where 7 = 3/2, S = —121/304, P = 3293/2299, /3 = —1181/2299, dj = 1 — (do + 2i) and do = —1/19. These constants 
are summarized in Table 1. In summary, the time and phase function up to IPN are 


t — tc = t^+tpN- (73) 

— 4>c = 4?jv + ^PN (^4) 


3 


tpN = — 


Sc^An 


= — 


+F(e,a2,/3,7:^)^4 
ISc^A"/" 


^J3(e,aj,/3,7,(5)Si 
— dTotAf, 


^PN = — 


304(G7Wc)^/3 

Sc^Ao 


F(e,d,/3,0,(5), 


leGMcV"^/^ 
+F{e,a2,^,0,S)B4 


J2F{e,a,J,0,S)B, 


-^ho 4 >Ar. 


(75) 


with Ao = (T(eo)/(27rt'o), To = PQ^^bp]s[{eo) and the function F(e, a,/3, 7 , 5) = Fi (^,/l, 7 , e“/a. 0 

The qualitative behavior of the orbital evolution is presented in Figures 1-3. 


^ It can be noticed that F{e,aN, (3 n^0,S) =2 for and ^p/v 
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TABLE I. Constants of the time and phase functions. 


6= 

N 

PN 

t — tc 

'7=2 
' 2 

^ = -ii 

^ = 3 
/ 2 

Q _ 1988 
^ 2299 

ao = a2 = 


0 

II II II 

7 = 0 

3 3293 

^ 2299 

ttO Yq ) ^1 Yq) ■ Yq 


eo 



FIG. 1. (color online). Comparison of the analytical and numerical results for the evolution of the eccentricity. The initial 
eccentricity is eo = 0.6, the masses of binary are rm = 10®(1 + z)Mq with redshift 2 = 1. The inspiral starts one year before 
the last stable orbit (LSO) calculated in the Newtonian order which corresponds to the initial frequency uq = 8.09/rHz. In 
the Newtonian case the dotted black line denotes the analytic, while the gray line the numeric solution. For the IPN orbital 
evolution the analytic and numeric solutions are denoted by the dotdashed and red lines, respectively. It can be seen that the 
perturbative solution is in perfect agreement with the numerical one. 


VIII. SUMMARY 


In our work we have investigated the orbital evolution and emitted radiation of binary systems on eccentric orbits 
up to IPN order. Both the time and frequency domain waveforms are presented in a simple form with use the 
generalized true anomaly parameterization. To express the time dependence of the waveforms the Hansen coefficients 
were generalized to IPN accuracy. Moreover, the radiation reaction problem and the evolution of the time and phase 
functions are given to IPN accuracy. 
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Appendix A: Tensor spherical harmonics 


Following the notation of [2^ the traceless, symmetric and unit basis tensors can be written as 

1 i 

ey® Cy) ± - (e^, 0 0 e^,), 

1 i 

= T 2 (Ca; 0 Cz + 0 e^;) - - (cy 0 0 Cy ), 

t -= ( Cg, 0 Ca; 0 Cy “t“ Z 0 ^2) * 

v6 


(Al) 
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f-f, (xIO^s) 



FIG. 2. (color online). The time function for various initial eccentricities. Newtonian and IPN expressions are denoted by 
solid (black) and dashed (red) lines, respectively. The masses of the components are rrii = 10® (1 + z)Mq with redshift z = 1. 
The inspiral time is set to 1 year before LSO calculated in the Newtonian order. 


(p-0^{x1OOOrad) 



FIG. 3. (color online). The phase functions for various initial eccentricities. Newtonian and IPN expressions are denoted by 
solid (black) and dashed (red) lines, respectively. The masses of the components are mi = 10® (1 + z)Mq with redshift z = 1. 
The inspiral time is set to 1 year before LSO calculated in the Newtonian order. 


The scalar harmonic tensors on this basis are given by 

l' 2 

Y , (A2) 


where {l'2m'm",lm) denotes the Clebsch-Gordan coefficients and T*"* is the conventional spherical harmonic. Then 
the electric and magnetic tensor harmonics can be expresssed as 


-^E2,lm 


i{i + i) 


2(21 + l)(21 + 3) 


T 


'2 l+2,lm I rri2 I 


Im 


(21 + 1) (2/ + 3) 




rj-iB2,lm _ _^ 


I (/ + !)(! + 2) 2 1-2,Im 

2 ( 21 - 1)(21 + 1 ) 


1-1 

21 + 1 


T 


2 




21 + 1 ' 


(A3) 

(A4) 


As an example we consider the tensor harmonics and T® 2,22 appearing in the Newtonian waveform. Using the 

relationship between the Descartes and spherical polar coordinates, 


Cx = Cr sin 9 cos + + eg cos 9 cos ip — e^ sin ip, 
Cy = Cr sin 9 sin ip + eg cos 9 sin ip e^, cos ip, 
Cz = Or COS 9 — eg sin 9, 


the tensor harmonics have the form 


rpE2,22 _ 


^B 2,22 


- Y^^[(l + cos^ 9) h+ + 2i cos01ix]e^*‘^, 


1 

16 V 27r 


[(3 + cos 20) hx — 41 cos01i+]e 


2iip 


(A5) 

(A 6 ) 

(A7) 
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where h+ = eg 0 eg — and hv = eg® etp-\- e^p® eg are the two independent polarizations. The tensor spherical 

harmonics up to 2PN are given in |18l |. 


Appendix B: Orbital parameters of the IPN dynamics 


The orbital parameters up to IPN order are given by 


with the quantities A — K are 


A 

B 

C 

D 

H 

I 

K 


{-A) 


3/2 


aR = 


et = 


B ’ 

B 

D 

A 

^ 2Co’ 

1 - 


( 1 - 

AD ^ 


1 1/2 


Co / 


et, 


eg = [1 + 

k = 


2BCo 
AD AI \ 

'mi) 


BCo 
3Gm 


ar(l ^ r ) 


et, 


2C A 3 E 

Gm(i + (7,-6)^). 

72 / E 

1 + 2(37,-!) — 

/ n o^ CMJ2 

(-317 + 8) 2 2 , 

- (l + (3??-l)-^ 

M V fj-c 

GMJ 

27,-4 - 

fiC‘‘ 

J 


ea^M'^ 


1/2- 


(5?, - 10) 


G^M^ 


(Bl) 


(B2) 


Appendix C: Hansen coefficients 

The Hansen coefficients are important functions of the celestial mechanics which are known for more than 100 years. 
The expansion of Hansen-coefficients is 

00 

exp(77770) = ^ exp(ikA4), (Cl) 

k— — co 

where r is the relative distance, a is the semimajor axis, (j) is the true anomaly and M is the mean anomaly using 
by the standard notations of celestial mechanics. The coefficients are called the Hansen-coefficients. Here the 

constants n and m are integers. The Fourier series representation of the Hansen coefficients is 

TT 

^n,m _— f Qxjp(im(j) — ikAi)dAi. (C2) 

ZTT 7 Va/ 
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The integrand can be transformation to other argument, specifically the eccentric and true anomalies with the use of 
the leading order Kepler-equation, 


dM , 

—— = 1 — e cos u, 
du 

(C3) 

dAi (1 - 

(C4) 

dv (l + ecos(/)) 


We introduce the complex variables x = expjc/i and y = expzu (z = expiA^ for the contour integral), then the 
relationship between the eccentric and true anomalies ([s^) with x and y complex variables is 


and one gets for the variable x 
and the mean anomaly Ai 


X — 1 1 + 13 y — \ 
x + l 1 — j3 y + V 

expti^ = 2/(1 - ^2/”^)(l - /3y)“\ 




exp iAi = y exp 
The integrand with eccentric anomaly u is given as 

(l + /3^)“”-' 




(C5) 

(C6) 

(C7) 

(C8) 


- 


27r 




— 1 \n+m+l 


ke 


du. 


(C9) 


(CIO) 


X (1-/32/)”-™+'exp y(j/-y') 

The integral can be extended to infinity as a series of the Bessel functions 

CXD 

A”’™ = (1 + /32) —1 ^ E-r,Uke). 

P— — 00 

The coefficients y”’"* for I ^ m and for I < rn can be expressed by the hypergeometric function F(a, 6; c; d) as 

— m + E 
I — m 

xF{l — n — 1,—m — n — l;l — m + 1;/3^). (Cll) 

The first description of this formula was given by Hill [2^. The other representation of the Hansen coefficients can 
be found in the work of Tisserand on celestial mechanics from 1889, 


j^n,m / Q\l — m ^ 3 ” 


A”’™ = 


'2s 


(l+^2)n+l 


s=0 


where 


Vs = 


Qs = 


Ps+k-m k^m { 

Ps /c < m J ’ 

Qs fc ^ m 1 

Q s+m — k k ^ Ul \ 


(C12) 

(C13) 

(CM) 


r—O 

s 

= E 


r—0 


n — TO + 1\ 1 / fcr' 
s — r ) r\ \2I3J ’ 

n + TO+l\l / kr\ 


s — r J rl \ 2/3 J 


(C15) 

(C16) 


and 
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Appendix D: Leading and half order waveforms 

The leading order waveform with the true anomaly (f> due to Einstein quadrupole formula are (the notation of ref. 
m for azimuthal polar angle is 7 = $) 


1 

1 — 




(Cm'*’ cos mcj) + sinTO(()) , 

m—0 

(Dl) 

5 

y^ cos m 4 > + sinm^) , 

m=0 

(D2) 


where 




= — — (1 + 3 cos 2$ — 2cos20sin^ $), 


= -(—1 + cos 20 — 5(1 + cos^ 0) cos 2$), 
C 2 ~^ = — (3 + cos 20) cos 2$, 
c^'*' = ~^(3 + cos 20 ) cos 2 <i), 

5e 

s^'*' = ~^(3 + cos 20 ) sin 2 <i), 

S 2 ~^ = — (3 + cos 20) sin 2<i), 

= — ^(3 + cos 20 ) sin 2 $, 


Cq^ = 2e^cos0sin2$, 


„Nx 


= 5ecos0sin2$, 


C 2 ^ = 4 COS 0 sin 2$, 
= e cos 0 sin 2$, 


Nx 


= —5 cos 0 cos 24), 


S2^ = — 4 COS 0 COS 24), 


„Nx 


= —ecos0cos24). 


(D3) 


The half order waveforms (denoted by the superscript H) are 

5 




hx W = 


16 


2^3/2 E (c)J) cos m(j) + sin m(j)) , 

(1 ~ 6 ) m =0 


m—0 

5 


(l-e2)3/^^ 


(c^ cos m(j) + sin m- 


4>) . 


(D4) 

(D5) 


where 


= —e (2 (—11 + 2e^ + (—1 + 6e^) cos 20) sin 0 sin 4) — 2e^(5 sin0 + sin 30) sin 34)) , 

= — - (2 (—44 — 19e^ + (—4 + 39e^) cos 20]) sin 0 sin 4? — 35e^(5sin0 + sin 30) sin 34)) , 

C 2 ^ = e((23sin0 — 5 sin 30) sin 4) + 15(5 sin 0 + sin 30) sin 34)), 

i (e^(—23 sin0 + 5 sin 30) sin 4) — 2 (I 8 + 7e^) (5 sin0 + sin 30) sin 34)) , 

= 5e(5 sin 0 + sin 30) sin 34), 

3 

cf+ =-e^(5sin0 + sin30)sin34>, 

i cos 4) ((42 — 72e^ + 175e^ cos 24)) sin 0 + (2 — 32e^ + 35e^ cos 24)) sin 30) , 

82 '^ = — e((23sin0 — 5 sin30) cos 4) + 15(5 sin 0 + sin 30) cos 34)), 

i (e^(23sin0 — 5 sin 30) cos 4) + 2 (I 8 + 7e^) (5 sin0 + sin 30) cos 34?) , 

= —5e(5 sin 0 + sin 30) cos 34?, 

= — ^e^(5 sin0 + sin 30) cos 34). 


(D 6 ) 
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Appendix E: IPN waveform 
1. The coefficients proportional to cos rruj) 

There are coefhcients proportional to cos m4> and sinm^ in the IPN waveform, Eq. (II3- Coefficients for cos 4^ and 
sin4(^ are (we have introduced the shorthand notations Cm<s> = costo$, 5^$ = sinmd), c^e = cosm0, SmO = sinm©, 
Ge = —5 + 4c2e + C 4 e and Xrj = — 1 + 3ry. For example, rnean plus polarization, proportional to sin3(/) and 

proportional to sin4(()) 

‘"4(c4) “ '^C4^GsXriTie , 

''4(s4) “ -^S44,GQXrirje , 

55 

''3(c4) “ ^^C4$GeA,,?7e , 

55 

'"3(s4) “ ^^■54$GeA^?7e , 

+ 125 ^ ^ 2 

''2(c4) “ ^g2'^4$GeA,,?7e , 

+ 125 ^ . 2 

*"2(84) “ ^g2'54$GeA^?7e , 

Ci(c 4 ) = ^C4$G'eA,,?7(764 + 135e^)e, 

^ 1 ( 84 ) = ^S4$GeA^77(764 + 135e^)e, 

= Y^C44.GeA,,?7(64 + 71e^), 

ciJ(84) = Y^S4$GeA^?7(64 + 71e^). (El) 

The coefficients proportional to cos 2(j) and sin 2(j) are 

^ 1 -j 

'^3(c2) = [291 + C4e - 377(27 + C4e) + 4(19 + 977 ) 020 ] rje'^, 

'^3(s2) = + '^40 - 377(27 + C 40 ) + 4(19 + 977 ) 020 ] 770 ^, 

^t{c2) = “^C2$(-983 + 114e^ + 37577 + 180^77 + C 4©(1 + 6 e^)A,, + 2o2©(-116 - 8 I 77 + 6e^(3 + r])))e^'q, 

^2{s2) = -^S2$(-983 + 114e^ + 37577 + 180^77 + 040(1 + 6 e^)A^ + 2o20(-116 - 8 I 77 + 6e^(3 + r])))e‘^T], 

o^^ 2 ) = i^*^ 2$(5948 + 99e^ - 3 OI 277 - 1089770 ^ - 040(4 + 81e^)A,, + 4 o 2 ©( 33 e^(-l + 77 ) + 4(76 + 81r])))er], 
c^g 2 ) = i^® 2$(5948 + 99e^ - 3 OI 277 - 1089776 ^ - 040(4 + 81e^)A^ + 4 o 20 ( 33 e^(-l + 77 ) + 4(76 + 8 l 77 )))e 77 , 
c+(^ 2 ) = ^C2<i.(508 + 0^(353 - 46677 ) - 22877 - O40(-4 + 19e2)A^ + 2020 (52 + 6 O 77 + e^(14 + 57tj)))tj, 

'^ 0 ( 82 ) = ^'S2$(508 + 0^(353 — 46677 ) — 22877 — C 40 (—4 + 19e^)A^ + 2c20(52 + 6 O 77 + e^(14 + b7r])))r]. (E2) 

The coefficients without (j) dependence are 

C 3 = [(5 - 1577)040 + 11(-11 + 77 ) + 4(29 + 77 ) 020 ] 770 ^, 

5 

cj = ^ [37 + 4(-9 + 77)020 - 777 + Xr^Cio] 0 ^ 77 , 

0 + = :^(1196 + 675e^ - I 9677 - 2650^77 + O40(-4 + 39e^)A^ + 4 o 20 (- 3 OO + 52r] + e^{-159 + 37ri)))er], 

C+ = ^?7(363 - SOe^ - 7377 - 100^77 + O 40 (-l + 6e2)A^ + O20(-364 - 8e2(-7 + 77 ) + 7677 )). (E3) 


18 


2. Coefficients proportional to sinm0 


The coefficients proportional to cos 4^ and sin 4^ are 


®4(c4) “ 22 ) 

®4(s4) “ 7^C4^GQ\r]e , 


32 

25 

®^(c4) “ 

25 

* 3 ( 84 ) “ -^^Ci<s>GQXririe , 

15 

®2(c4) Y 2 ’ 

15 

^2(s4) 22 ^4:^GQXrjT]€ , 

s^c 4) = Y^GeS44.A,,(52 + 9e^)?7e, 


«i( 84 ) = -Y^GeC4<i.A^(52 + 9e‘^)r]e. 


(E4) 


The coefficients proportional to cos 2^ and sin 2^ (coefficients of sinmt/)) are 


'®3(c2) = ^s 2$(-293 - 92c2e + C4e + 37?(29 + 4c2e - C4e))r7e^, 

*3(82) = “^C2$(-293 - 92c2e + C4e + 3?7(29 + 4c2e - C4e))?7e^, 

^2(c2) ~ Yg'®24>(~259 + 8177 + C 4 o(—3 + 2e^)A^ + 2e^(19 + 377 ) + 4 c 2 ©( 2 (—10 + 77 ) + e^(3 + , 

^^( 82 ) “ —16^^*^^— 8 I 77 + C 40 (—3 + 2e^)A^ + 2e^(19 + 377 ) + 4 c 20 ( 2 (—10 + 77 ) + e^(3 + Ti)))rje^, 


^i(c 2 ) — 64 ^^^^—20477 + 5c 4©(—4 + 3e )Xri H- lie (1 + 2 I 77 ) H- 4c20(—4(17 H- 77 ) 4- e (1 + 1977)))77e, 

4 ( 82 ) = -^S2$(-900 + 20477 + 5c40(-4 + 3e2)A^ + lle"(l + 2 I 77 ) + 4c20(-4(17 + 77 ) + e2(l + 1977)))77e. (E5) 

The coefficients proportional to cos 4^ and sin 4^ are (we have introduced the notation Hq = cq — C 30 ) 


^4(c4) ^'^44>-^0'^?7^e , 


4(c 4) = -7S4$-ff0A^77e‘‘ 


55 


^3(c4) ~ ^S44)iJ0A,,77e , 

55 

* 3 ( 04 ) “ —-^Ci^HQXrfTie , 

125 

‘'2(c4) “ ^^S44>-ff0A,)77e , 

125 ^ 2 

= -^^C4$i70A^77e , 
1 


^2(c4) 


4(c 4) = Y^S44.i?0A^77(764 + 135e )e, 
■^ 44 ) = -Tc 4 <i.i 70 A^ 77 ( 764 + 135e^)e. 


(E 6 ) 



The coefficients proportional to cos 20 and sin 20 are 

'^3(c2) = [(3 - 9?7)c3e + 11(-17 + 3?7)ce] 

*3(02) = [(3 - 9?7)c3e + 11(-17 + 3?7)ce] rje^, 

'^ 2 (c 2 ) = [(-17 + 6e^)A^C3e + (-625 + 1597? + 6e^(13 + r]))ce] r/e^, 

'® 2 (c 2 ) = [(-17 + 6 e^)A^C 3 e + (-62 5 + 1597? + 6e^(13 + r 7 ))ce] rye^, 

^i(c 2 ) = 4^2$ [5c3e(-28 + 9e2)A^ + (3(7 + 15577)e2 + 4(-931 + 321ry))ce] rye, 
^i(c 2 ) = "4^2$ [5c3e(-28 + 9e")A^ + (3(7 + 155ry)e2 + 4(-931 + 32l7y))ce] rye. 
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